/*==========================================================================================    * Date: Sept 2018
	* Purpose: To construct the Distance Matrix of the MOLA points that merged with residential points from the NAPP Census 
	* DATASETS Used:	
		- "ll_point1881_soc_civ_near_napp.shp" comes from merging in ArcMap of "MOLA\ll_points1881_soc_civ_near.shp" and "MOLA\MolaModified_NAPP_ALL_JAG1new.txt" (created by MolaNAPP16v3.do)
	* DATASETS Created:
		- "ll_point1881_soc_civ_near_napp.dta"
		- "points1881DistMatrix.dta"
==========================================================================================*/
	
clear all
set maxvar 20000
set more off

gl MOLA "/Users/`c(username)'/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/MOLA"

// Transform shapefile into STATA

shp2dta using "$MOLA/LL_POINTS1881_Soc_Civ_Near_NAPP.shp", database(ll_point1881_soc_civ_near_napp) coordinates(londoncoord2napp) genid(id) replace


use ll_point1881_soc_civ_near_napp, clear

keep OBJECTID_1 POINT_X POINT_Y
rename POINT_X longitude 
rename POINT_Y latitude 

levelsof OBJECTID_1, local(id_point)

foreach id1 of local id_point {
	cap drop d`id1'
	gen d`id1'=0
}

local id_left: list id_point | id_point

foreach id1 of local id_point {
	di "ID1 : `id1'"
	qui sum latitude if OBJECTID_1==`id1' 
	qui local lat`id1'=r(mean)
	qui sum longitude if OBJECTID_1==`id1' 
	qui local lon`id1'=r(mean)
    qui local id_left: list id_left - id1
	
	foreach id2 of local id_left {
		qui sum latitude if OBJECTID_1==`id2' 
		qui local lat`id2'=r(mean)
		qui sum longitude if OBJECTID_1==`id2' 
		qui local lon`id2'=r(mean)
		qui local dist=sqrt((`lat`id1''-`lat`id2'')^2+(`lon`id1''-`lon`id2'')^2)
		qui replace d`id2'=`dist' if OBJECTID_1==`id1'
		qui replace d`id1'=`dist' if OBJECTID_1==`id2'
	}
}

keep OBJECTID_1 d*

save "$MOLA\points1881DistMatrix.dta", replace
